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We study an instantaneous bimolecular chemical reaction in a two-dimensional chaotic, incom- 
pressible and closed Navier-Stokes flow. Areas of well mixed reactants are initially separated by 
infinite gradients. We focus on the initial regime, characterized by a well-defined one-dimensional 
contact line between the reactants. The amount of reactant consumed is given by the diffusive flux 
along this line, and hence relates directly to its length and to the gradients along it. We show 
both theoretically and numerically that the probability distribution of the modulus of the gradient 
of the reactants along this contact line multiplied by k2 does not depend on the diffusion k and 
can be inferred, after a few turnover times, from the joint distribution of the finite time Lyapunov 
exponent A and the frequency i. The equivalent time r measures the stretching time scale of a 
Lagrangian parcel in the recent past, while i measures it on the whole chaotic orbit. At smaller 
times, we predict the shape of this gradient distribution taking into account the initial random 
orientation between the contact line and the stretching direction. We also show that the probability 
distribution of the reactants is proportional to k 2 and to the product of the ensemble mean contact 
line length with the ensemble mean of the inverse of the gradient along it. Besides contributing 
to the understanding of fast chemistry in chaotic fiows, the present study based on a Lagrangian 
stretching theory approach provides results that pave the way to the development of accurate sub- 
grid parametrizations in models with insufficient resolution for capturing the length scales relevant 
to chemical processes, for example in Climate-Chemsitry Models. 



I. INTRODUCTION 

Chemical reactions in the stratosphere have been 
shown to be sensitive to the numerical spatial resolution 
when the chemistry is fast compared to advective pro- 
cesses (P, 111)- It was proposed by [l| that the product 
concentration of the deactivation of polar vortex chlorine 
by low latitudes nitrogen oxide at the edge of the strato- 
spheric Northern hemisphere winter time polar vortex 
scales like k^*^*) , with being k the reactant diffusion. Later 
on, [3] argued that p{t) can be expressed as 1 — D{t)/2, 
where D{t) is the box counting fractal dimension of the 
contact line between the reactants. Here we focus on 
the initial regime of a instantaneous bimolecular chemical 
reaction in a two dimensional Navier-Stokes flow charac- 
terized by chaotic trajectories (this provides an idealized 
framework for isentropic dynamics in the stratosphere). 
By definition, the initial regime is characterized by a well- 
defined one-dimensional contact line (i.e. D = 1). The 
reactants are initially separated by infinite gradients. 

In a previous work ([j]) dealing with this regime, we 
have shown that the ensemble mean reactant concentra- 
tion time derivative scales like k2 and can be predicted 
accurately from the Lagrangian stretching properties of 
the flow. Here we investigate the statistical properties of 
the chemical production and of the reactants concentra- 
tions. 

In section II, we explain how the study of an infinitely 
fast chemical reaction A + B — > C simplifies into the 
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study of a passive tracer whose concentration field </> is 
defined as the difference between the concentrations fields 
of the two reactants A and B . The rate at which the 
reactants disappear is the diffusive flux of (j) along the 
contact line, and hence depends on both its length and 
the gradients along it. In section III, we give some the- 
oretical relations between, on one hand, the contact line 
length and the gradients of the reactants along it, and, 
on the other hand, the Lagrangian stretching properties 
of the flow. In section IV, we focus on the probability 
distribution of the reactants concentration. Section V 
describes the Lagrangian stretching properties of a two- 
dimensional Navier-Stokes flow and section VI presents 
some numerical simulations to test the theoretical results 
of sections III and IV in the flow introduced in section 
V. 



II. FINITE TIME LYAPUNOV EXPONENTS 
AND CHEMICAL PRODUCTION IN A 
CHAOTIC FLOW 

Wc consider the bimolecular chemical reaction A + 
B — > C in stoechiometric quantities. One molecule of A 
reacts locally with one molecule of B to give one molecule 
of C. As a result, the field (p = Ca — Cb, defined as the 
difference between the reactants' concentration fields Ca 
and C_B , is independent of the chemistry. If, in addition, 
the reaction is instantaneous, A and B cannot coexist at 
the same location. Consequently, the fields Ca and Cb, 
and their spatial average over a closed domain denoted 
by an overbar, can be retrieved from (j) as follows (0, 0]): 
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If A and B are separated by a contact line C = {x|0(x) = 
0} of dimension one oriented in a counterclockwise di- 
rection such that it encloses reactant A (domain Da), 
the time derivative of the reactants in an incompressible 
closed flow is: 
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where A is the total area of the domain and n a unit 
vector normal to £{t) pointing outward from Da- We 

call — the chemical speed. Furthermore, on every 
line element dl of C{t), the quantities f V0 • ndl of A 
and B arc consumed. As a consequence, knowing the 
length of C{t) and the probability distribution of IV0I 
along it, gives a comprehensive statistical description of 
the chemistry in the domain. 

The equation for the passive tracer (jj is: 



^ + u • = kV^(; 
ot 
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where u is the flow. If the trajectories are chaotic, equa- 
tion ([2]) allows to link the chemical speed to the La- 
grangian stretching properties of the trajectories as cap- 
tured by the finite time Lyapunov exponents (FTLE), 
defined as the rate of exponential increase of the dis- 
tance between the trajectories of two fluid parcels that 
are initially infinitely close. If Sl{t) is the distance be- 
tween two parcels that start at x and x -I- ^Iq, then the 
FTLE X{x,t) at x over the time interval t is 



M^,t) = Tniaxj In— — I , 
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where the maximum is calculated over all the possible 
orientations a of SIq. The unit vector with the orienta- 
tion ■0_|_(x,t) of SIq at the maximum defines a "singular 
vector" i/j_|_(x, i) = (cos i/'-i- , sin 7/1+ ) . In the inviscid limit, 
it can be shown, with the conservation of tracer concen- 
tration for Lagrangian parcels, that A(x, t) is also the 
rate of exponential increase of a gradient initially aligned 
with the unit vector 'i/j_(x, t) = {— sinip-^., cos ip+) per- 
pendicular to ^p-\-■ 

We can calculate FTLE and singular vectors in an 
incompressible flow using the velocity gradient tensor 
S = Vu(X,t) along a trajectory X(x, <). The distance 
51 between two trajectories initially infinitely close is 
solution of ^ — S{t).6l = and is given by 61 = 
'M.6l{t — 0) where the resolvent matrix M is solution 
of ^ - S(i).M = 0. The finite time Lyapunov expo- 
nent A(t) is given by the log of the largest eigenvalue of 
[M-^M]2t , with the associated eigenvector. 



In ergodic chaotic dynamical systems, it has been 
shown that the FTLE converge to an infinite time Lya- 
punov exponent independent of the initial position x, 
while the singular vectors converge to the forward Lya- 
punov vector = (cos sin ^i^) that depend on x 
(Osedelec theorem, Q). The convergence of the Lya- 
punov exponent is very slow and algebraic in time ([6|) 
while the convergence of the Lyapunov vector is much 
faster, and typically exponential (Q). A discussion spe- 
cific to high Reynolds number two-dimensional Navier- 
Stokes flows is available in [S]. Here we will only take 
into account the time dependence of the FTLE, assum- 
ing the singular vector field is only a function of space, 
not of time, equal to the forward Lyapunov vector field. 
This assumption will allow us to link the evolution of 
a gradient along a trajectory to the Lagrangian strain- 
ing properties of the flow, while taking into account the 
diffusion (see (fTO|) below) 

An element SIq — \dlo \ (cos a, sin a) of the contact line 
at the initial time is advected at time t into an element 
dl whose squared norm is 



\dl\^ =dllM'^MSlo 
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cos^(?/)+ — a) + e sin (^/j. 



a)] 
(5) 

Noting that the angle a is random, we can show, aver- 
aging over A, tp+ and a, that the total ensemble average 
length of the contact line is (brackets are for an ensemble 
average) : 



lo ir ^dlPx {t, 0Ve2"cos2 7 



e 2'* sin^ 7, 
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with Lq the initial length of £ and P\{t,l) the time de- 
pendent probability density of A. Henceforth, the inte- 
gration bounds over I and 7 will be implicit and the same 
as in ([6]). This expression is valid when the diffusion is 
not taken into account. For large times, when two fila- 
ments are brought together at a distance smaller that the 
diffusive cutoff, they merge under the action of diffusion. 
We think the time span of this regime to be well approxi- 
mated by the mix-down time scale Tmix ~ ^ In. from 
the largest scale of the flow to the diffusive cutoff ( 0] ) . 
If the strain S is an estimation of A, which is the case for 
the flow we study section V and VI (flgure[T|), and if ^ 
is an estimate for the integral time scale T of the flow, 
which is also true in our flow, we get that 
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where Pe is the Peclct number, Re the Reynolds number 
and Pr = ^ the Prandtl number. Our previous work 
3 has shown that ^ is very accurate on time scales of 
the order of Tmix in the flow we describe in section V. 

Assuming the stationarity of the singular vectors, as 
explained previously, and noting the existence of a local 
1-D solution Xt(-2;), as in [loj . in the direction z normal 
to the contact line, of the advection-diffusion equation 
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([3]) written in the co-moving frame with a contact line 
element, Q showed that, for an initial gradient profile 
of tracer AqS{z), with S the Dirac delta function, this 
solution is: 



Xt{z) = 
with G 



^oErf(^ 



e2^t cos2(iI>^— Q)+e-2At sin2('I'^— q) 
re^^* cos2(iI>^— a)-|-Tsin^('I'+— a) ' 



(8) 



The concentration Aq is the initial concentration of reac- 
tants A and B in their respective domain. The function 
Erf is the gauss error function defined on R as follows: 
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dt. The two quantities t and f arc: 
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The time r has been introduced through the wavenum- 
ber growth along Lagrangian trajectories by [ll| and was 
called an equivalent time by [l^ . Because the trajectories 
are chaotic, ^ is the stretching rate in the recent past (i.e 
approximately over the last correlation time). [l2| have 
argued that when the correlation time of the stretching 
is much smaller than the time scale of the decay of the 
mean Lyapunov exponent, t becomes independent of A 
at large times and that its probability distribution con- 
verges to a time independent form. The time r is also an 
equivalent time that measures the stretching rate in the 
early part of the trajectory. As a consequence, we expect 
T and T to have the same statistics, to be asymptotically 
equivalent as t — )■ (typically for times smaller than the 
correlation time of the stretching) and to become inde- 
pendent at larger times. From ([9]), we can see that the 
gradient along the contact line |V0£| = ^^^-^| 



is: 



(10) 



For times such that t ^ j^, the sine terms in (|TT|) can be 
neglected: the gradients are equilibrating with the flow 
and the time asymptotic form of G is (see ([8])). As 

a consequence, for t ::S> -jy w ^, PcxitTg) is asymptoti- 
cally equivalent to Pg,c,<x>'- 



PG,£(t,g) -i"G,£,oo(i,.9) 
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JdlPx{t,l)e^' ' ^^^^ 



where P ^ \ is the joint pdf of and A. It is worth not- 
ing that if r and A were independent, which is expected 
at the very long times, the pdf of g along the contact 
line would be equal to the pdf of This can be seen 
writing the bivariate density P i as the product of its 



marginal densities in ([131 



IV. PROBABILITY DISTRIBUTION OF THE 
REACTANTS CONCENTRATIONS 

The reactants are in stoechiometric quantity and their 
initial concentration in their respective domain is Aq. As 
a consequence, it follows from ([T]) that the pdfs of Ca, 
Cb and |0| are the same. The corresponding random 
variable will be noted $. Our objective is to understand 
its dependence with k and with time, as well as its shape. 

The profil of the tracer gradient close to the contact 
line is given by (|5]), expression that should be valid as long 
the contact line is well-defined with a curvature larger 
than which is expected as long as t < Tmix- We 
assume that an e ^ 1 can be chosen such that, for all 
members, the length 



III. PROBABILITY DISTRIBUTION OF THE 
REACTANT GRADIENTS ON THE CONTACT 
LINE 

The distribution of | V(/)£ | -^^^ can be inferred form the 
distribution of G (eq. [TO)) and A and does not depend 
on K. Its probability density function (pdf) along the 
contact line Pg.c is given by: 
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where we have introduced the joint pdf Pg,\ of G and 
A. Considering an orbit characterized by (A, r, r, 4'-)., a), 
expression pT|) can be derived noting that ^^^^ |V(/)£| is 
equal to G(t, A, r, r, 4'+, a) on a fraction of the contact 
line 
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2xrH^o(l-6)) = ^Erf-i(l-e) (14) 



> (5t > The range [-^ ^ 



This 



satisfies La = 

is where \xt\ takes values smaller than Aq{\ — e) 
can be seen noting that xt is a monotonic odd increas- 
ing function with xti,'^^) = ±^o(l ^ e)- In fact, we 



have 



<C 1 because 



'Pe is very large by 



assumption (the Peclet number in the simulations pre- 
sented in section VI will be of the order of 10* to 10^). 
However, ^ is in general different form S: the former 
measures a Lagrangian stretching rate in the recent past 
while the latter measures an Eulerian stretching rate (for 
a comparison in the flow considered in section VI, one 
can refer to figure 1, noting that the strain is the Lya- 
punov exponent at small times). It is possible to choose 
5t far from both La and ^/kt because we are considering 
a well-defined contact line between areas of well mixed 
reactants. As mentioned earlier, this is expected to be 
true for t < Tmix- 

If we consider a profile Xt of tracer around the con- 
tact line in the range [— ^, the cumulative distribu- 
tion function (cdf) of $, i.e. the probability of having $ 
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smaller than a given value 4>, is 



xt\4>) 
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'xr'(^o(l-e) Erri(l-e) 

(15) 

and stands for values of </> in the range [0,Ao(l — e)]. 
We multiply by ^"^'^^"^^ to obtain the cdf F^f, in the whole 
domain of area A because values in the range [0, Ao(l — e)] 
are achieved on a sub-domain of area {6t){L): 
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Using equation (|T4|) to calculate (St), we have: 



(16) 



i^*(0) = ^^Wi^) Erri(^) for G [0,^o(l - e)]. 

(17) 

Finally, the probability density is the derivative of the 
cdf with respect to 0: 

(18) 

A direct consequence of well mixed reactants away from 
the well-defined contact line is the term Erf~^ (deriva- 
tive of the inverse of the Gauss error function). It ex- 
presses that the shape of the pdf is an increasing func- 
tion and can be directly related to the profil of the tracer 
field close to the contact line. The denisty P$ is pro- 
portional to \fti{L) because the area where the field 
takes non-trivial values (i.e significantly different from 
the initial value ^o) is proportional to ■\/k(L): its length 
is (L) while its width is controlled by diffusive processes. 
Finally the term (^) depicts the effect of the mean gradi- 
ent, with a decrease of the gradient along the contact line 
explaining an increase in the probability of small values 

of !</>!. 



V. STATISTICS OF THE LAGRANGIAN 
STRETCHING PROPERTIES IN A 
TWO-DIMENSIONAL NAVIER-STOKES FLOW 

The numerical model integrates the following vorticity 
equation using the pseudo-spectral method: 



dijj 
'dt 



u • Vuj = F - Rqoj + vV uj. 



(19) 



where w is the vorticity, F a forcing term at wavenum- 
ber 3, i?o the Rayleigh friction and v the viscosity. 
The equation is integrated in a doubly periodic box 
[— TT, tt] X [— TT, tt] on a 512 x 512 grid. The integral time 

scale of the flow T = 



2 



■J, where the brackets 



stand for an ensemble average, will be used to normalize 
the time axis. The Reynolds number is of the order of 
10". 
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FIG. 1. Probability density of the Lyapunov exponents (top 
left), of the inverse of the equivalent time t defined in @ 
(top right) for < f < 25r. We note that the density of A at 
t = 0.25T is roughly the density of the strain. At the bottom, 
their joint density at t = 4r (left) and t = 25T (right). 



Trajectories are computed using a fourth order Runge- 
Kutta scheme with a trilinear interpolation on the veloc- 
ity field. On each trajectory, we integrate the resolvent 
matrix M such that ^ = SM with S = the veloc- 
ity gradient tensor along the trajectory and M(t — 0) the 
identity matrix. The largest eigenvector of the symmetric 
positive matrix *MM is e^"^*, where A is the Lyapunov 
exponent on the trajectory at the finite time t. This 
method is described in more detail in [l^. We also cal- 
culate the equivalent times r and r through a numerical 
integration of The trajectories are computed for hun- 
dred realizations of the flow, each realization spanning 25 
turnover times. This gives access to the statistics of A, 
T, T and G involved in equations ([TT]) and ([T^ . 

Figure [1] shows the time evolution of the pdf of A and 
i . The pdf of A converges to the pdf of the strain Pg as 
t — > because the strain is the FTLE on each chaotic 
orbit for an infinitely small time (see (|3])). The pdf does 
not evolve much during the first turnover time, as the 
correlation time is expected to be of the order of T, or 
larger. Then, the variance of the FTLE decreases while 
the density shifts toward smaller values. The peak of the 
density saturates at Xmax ~ 0.02, which we think is a 
rough estimate of the infinite time Lyapunov exponent. 
For a more detailed description of the FTLE, the reader 
can refer to [l3l in chaotic flow, to [1] in two-dimensional 
turbulence and to [ll, [H, [l^ in geophysical flows. For 
times smaller than one turnover time, the pdf of ^ shifts 
toward smaller value, its shape being only slightly af- 
fected. This can be interpreted assuming that A does 
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not evolve much on trajectories and can be estimated 
by the strain S where the trajectory originates. Hence, 
T w 25/(1 - e-2St)^ which gives 7 w A + 7 for < < ^. 
The pdf P{t,x) of i is then approximated by Ps{x — j) 
where Ps is the pdf of the strain. In addition, as we 
expected in section II, we observe that the density of ^ 
converges to a time independent form. 

Figure [T] also shows the joint pdf of A and ^ ai t = 4T, 
well within the time range where we think (|lip should 
be a satisfying description of the gradient pdf. The joint 
pdf at i = 25T shows this dependence is still impor- 
tant at large times. Previous studies (e.g [Tl|) have as- 



sumed the independence between A and r at long times. 
This is relevant in simple chaotic flows. However, two di- 
mensional Navier-Stokes flow exhibit coherent structures 
(vorticies, filaments of vorticity, etc..) probably making 
the Lagrangian correlation time dependent on the trajec- 
tories. In particular, very long correlation times could be 
associated with trajectories trapped in vorticies, where 
the stretching rate is particularly weak. This could ex- 
plain the strong dependence between large values of r 
and small values of A. 



since the points are not equidistant. The probability den- 
sities obtained after normalization of these histograms 
are noted PcPr- 

On figurc[3]wc plot PcPr for different Prandtl numbers 
and different times. We observe: 

• When time becomes shorter, as observed from t = 
4r to t = 0.25T, the independence of Ge.Pr to Pr, 
as predicted in sections HI and IV, is not verified. 
As a matter of fact, the gradient cannot be consid- 
ered as infinite at the initial time because of the 
finite grid size of the numerical model. This effect 
can be quantified: we can solve the advection dif- 
fusion equation in a Lagrangian co-moving frame 
with a contact line element for the tracer profile 
around this line with the initial condition on the 

^° " ^ {5q is a length cor- 



gradients %L=o 



dz lt=o 2(5oV^ 

responding to a grid point). We get 



with G^- = 



e^^* cos^(iI'+— sin^ (>!'+ — a) 



(21) 
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VI. NUMERICAL RESULTS 

A. Gradients along the contact line 

The numerical simulations are performed for eight dif- 
ferent Prandtl numbers Pr = = 2' for < i < 7. For 
each one wc run an ensemble of 34 simulations integrating 
equations (O and (|T9)l in the periodic box. Each member 
is defined by its initial condition on the flow, taken as the 
vorticity field every turnover time of a long time simula- 
tion of the statistically stationary flow solution of ([T^ . 
For each member, we use the following initial condition 
on the tracer: 

(j){x, y,t^O)^ 2Ao{H{x) - i) for {x, y) G [-tt, Tr]^, 

(20) 

where H is the Heaviside step function. In other words, 
in one half of the box, Ca = ^0 and Cb = 0, and in the 
other half Ca = 0, and Cb = Aq- A and B are separated 
by initially infinite gradients. This is not exactly true in 
the numerical integrations because of the finite resolution 
of the model and has to be kept in mind for an accurate 
interpretation of the numerical results. 

For each member we determine the coordinates of the 
contact line L with a time increment ^ using the library 
DISLIN ([l3)- We calculate the modulus of the gradient 
of (j) at each of these coordinates using bilinear interpo- 
lation. We then multiply it by -^^^ in order to obtain a 
physical quantity that scales like G in equation ([8]) and 
which we name G^^pr- A weighted histogram of G^^pr 
is calculated every time increment ^ using as weight the 
length enclosed by three consecutive points of C centered 
in the point where Gg.Pr is estimated. This is necessary 



where G^ now depends on diffusion through the 

term (this expression has to be compared to G 
in ([HJ). The initial gradient cannot be assumed in- 
finite when the diffusive cutoff y^ltr is of the order 
of the grid size So. The time scale Tpr where this ef- 
fect is important can be evaluated by comparing the 



two terms and [re^"^* cos^ (^1^+ — a) + r sin^ (^f^ — 
a)] in the denominator of G^. Approximating A 

with S, we find Tpr ~ ln(l + -^2.Pr), which gives 

respectively Tpr/T = 0.12; 0. 25'; 0.45; 0.75; 1.2; 1.6 
for Pr = 4; 8; 16; 32; 64; 128. This calculation is 
consistent with the numerical results presented on 
figure [21 the pdf P^Pr being independent of Pr for 

t>Tpr 

u At t = 4T, the densities PG,Pr=i2s a-nd PG,Pr=64: 
are different from PG.Pr<32- Because of the fi- 
nite resolution of the model, gradients at the high 
end of these distributions cannot be resolved. This 
explains the strong asymmetry of these densities. 
This effect is also observed at other times: in par- 
ticular, it explains the kick at the very high end of 
the densities Gpr>32 at t = 0.25T. 

• At t ~ IT, the densities PcPr are not anymore 
independent of Pr because we are at t > Tmixi^) 
for all the Prandtl numbers 

In order to remain consistent with the infinite initial 
gradient hypothesis, we will only consider small enough 
Prandtl numbers. We compare on figure |3] the numerical 
results to the theoretical predictions Pg.c ([III): Pg.c.oo 
(US]), and to the pdf of 
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FIG. 2. Probability density PG.i<Pr<i28 of Ge,Pr = a^\''^4'c\ where |V(^£| is the modulus of the gradient of <l> along the 
line C = {x|(/)(x) = 0}. These pdf are obtained from an ensemble of 34 direct numerical simulations and plotted here for 
t = \T, \T, T, 4T and IT 



• For t < T, only Pg.c shows some success in pre- 
dicting PcPr- This is consistent with (jlip which 
states that the effect of 7 cannot be neglected at 
times smaller than T. 

• For t = 4T, Pg.c a-nd Pg.c.oo sue much more sim- 
ilar because the contact line elements have equili- 
brated with the flow (their orientation does not de- 
pend anymore on their initial orientation a). The 
agreement between Pg.c and Pg.Pt is very good. 

• For t = 7T, Pg,c and Pgx.oo are even closer. 
Pg,c,oo fails to predict PG,Pr=i but performs rea- 
sonably for PG,Pr=32- Actually, an estimate of 
Tmico from gives 4.5T for Pr = 1 and 6.5T for 
Pr = 32, which is consistent with the discrepancy 



between PG,Pr=i and PG,Pr=32 at t = 7T. 

It is worth observing that the pdf of fails to predict 
Pg.Pt because of the dependence of A with r, which is sig- 
nificant even at very long times (section Vj3- We observe 
that the time scale for the much simpler Pg,c,oo to be- 
come a good prediction for the gradient, i.e the timescale 
for Pg.c to converge to Pg,c,oo, seems to be of the order 
of Tmix ■ This non-trivial behavior may be determined by 
the dependence of t with r ([HI). If wc had t ^ r, which 



^The similarity between the pdf of and PQ p^—i{t = IT) (figure 
Isjl is not explained by our theory, as showed by the significant 
difference between the pdf of -j^ and Pg,C- 
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FIG. 3. Comparison between Pc.Pr obtained from the direct numerical simulations and the theoretical predictions Pc.c 
(eq. Ilip . Pg,c,oo (eq. I13p and the pdf of (t is defined in ((Ojl), obtained from the calculation of the Lagrangian stretching 
properties of the flow calculated with the trajectories. We have only plotted the curves PcPr corresponding to direct numerical 
simulations consistent with the infinite initial gradient hypothesis. 



is expected as i — > 0, this convergence would have been 
of the order of ~ 5^ ~ ^ (equations (fTO)) and (fTTj) '). 
It is one order of magnitude longer. The reason could lie 
in the fact that r and t quickly become independent. 

Equation dTTI) . where the joint pdf of A and G is re- 



placed by the joint pdf of A and G^, is a fair approxima- 
tion to Pc.Pr <32 at any time (not shown). This extends 
our results to finite initial gradients, but concentrated 
on length scales not large compared to the diffusive cut- 
off of the flow, such that a Lagrangian straining theory 
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FIG. 4. Numerically determined pdf Qpr of = Ao — \(f>\ 
at different times t = \T, T, 2T, AT and IT and for Prandtl 
numbers Pr = 1, 4, 8 and 32. Log-log scale. 



FIG. 5. Numerically determined pdf Qpr oi 4> = Ao — |<^| 
multiplied by v^P^ for Pr = 1, 4, 8 and 32 and f = iT, T, t = 
4r and IT. Log-log scale. 



approach remains possible. 



B. Reactants' fields 

We calculate the pdf Q Pr oi (p = Aq ~ using our 34 
simulations ensemble for the entire range of Prandtl num- 
bers. On figure m we show Qpr for Pr — 1,4,8 and 32 
and for 0.25T < t < IT. Our objective is to compare the 
theoretical prediction to these numerical results. 

Figure [S] shows V PrQ Pr for a wide range of Prandtl 
numbers over the first few turnover times. The depen- 
dence in ^/K predicted in (fT8|) is well achieved, except 
at small times {t = 0.25T) especially for small diffusion 
(Pr = 32), i.e. when the infinite gradient assumption is 
again violated in the numerical simulations. It also fails 
&tt = 7T because T > T^i^. 

Next, we compare the time dependence predicted in 
([T5| to the numerical results. For Pr = 8, we show the 
product of Q Pr with -s/Pr-^^^^ on figure [6] The con- 
tact line length (L) is calculated form P\ using ([6]), and 
(1/G) from the integral ^Pc.cit, g)dg where Pc,c is 
defined in pT|) . The densities P\ and Pg,c are deter- 
mined as explained in section V from the computation of 
Lagrangian trajectories. The curves converges together, 
except at very small times for values of |0| close to (0 
close to 0). As expected, ((T5)) does not work either for 
t > Trnix, as shown by the curves t — 7T and t = 12T. 

The shape of the reactants' pdf shown as 
Qpr^sVP^^ on fi gure [5] is very well reproduced 
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1 1 




t=0.25T • : 




t=T 




t=2T 




t=4T 




t=7T 




t=12T 




\. Theoretical prediction : 


1 L 





1e-04 0.001 0.01 

Ao-|(^l 

FIG. 6. Time evolution of Qpr^s^/Pr^^^ . (L) and 
G are defined respectively in equations ((5]) and ([S} and 
are calculated from trajectories, as described in sec- 
tion V. The red curve (theoretical prediction) corresponds 

Erf^^ (^^A^)' 'where Erf is the Gauss error function. 
Log- log scale. 



9 



'XJ^^^^ ^ ( "^Ao ) ' '^t^i'^t^ confirms the theoretical 
prediction (|T8|) . 

VII. CONCLUSION AND DISCUSSION 

We have considered the early regime of an instanta- 
neous chemical reaction in a two-dimensional Navier- 
Stokes flow, for segregated reactants initially separated 
by infinite gradients. The time scales considered here arc 
shorter than the mix-down time scale from the integral 
length scale to the diffusive cutoff {Tmix given in (O). As- 
suming that the singular vector associated with a FTLE 
does not depend on time and is equal to the forward 
Lyapunov vector, we have adopted a Lagrangian strain- 
ing theory approach and showed numerically its success 
in predicting (a) the statistics of the diffusive flux of fast 
reacting chemicals along their interface (fTTj) . and thus the 
chemical speed, (b) the probability distribution of the re- 
actants ([T5)) . We have put the emphasis on the effect of 
the reactants diffusion, showing (a) that the distribution 
of the gradients along the contact line, rescaled by k2, 
does not depend on k (b) that the distribution of the 
reactant concentration, when it is not too close to is 
proportional to k^. 

At the very early stage, for about one turnover time, 
predicting these statistics requires to know the joint 
statistics of G(A, t, r, 7 = — a) (equation [5]) and A, 
which is given by the joint statistics of A, t and r (7 
is a random angle because a is a random angle, and is 
thus independent of the other variables). At moderate 
time scales (a few turnover times), the knowledge of the 
joint pdf of A and t is sufficient. More work is needed to 
understand these distributions, i.e the marginal distribu- 
tions of A, T and t and how these three variables depend 
on each other. 

Previous studies ([IB, [13) investigating the two- 
dimensional mixing in the upper layer of the ocean sug- 
gest that it might be possible to recover the distribu- 
tion of A from the distribution of the strain (which is 
A as i — >■ 0), a more accessible Eulerian quantity, and 
from the time evolution of the mean Lyapunov expo- 
nent (A) . Specifically, they observed that the distribution 
of -TTv does not change significantly in time over several 



turnover times and is given by a Weibull distribution. We 
have observed the same behavior in our fiow (not shown), 
with the distribution of the strain being very well ap- 
proximated by a Weibull distribution of shape parameter 
1.^. The time evolution of (A)_in a chaotic ergodic fiow 
was theoretically predicted by [a] , a prediction which was 
shown to perform well in the mixing layer of the ocean 
f[l3j). Nevertheless, the mechanisms in two-dimensional 
Navier-Stokes flows (or in similar geophysical flows) be- 
hind the evolution of the finite time Lyapunov exponent 
distribution toward smaller values (figure [IJ were not ad- 
dressed in the literature to our knowledge. Although this 
aspect deserves more investigation, we can speculate that 
parcels tend to stay longer in area of low strain than in 
areas of high strain. For example, a parcel can be cap- 
tured in a vortex for a long time but travels very quickly 
in areas of high strain. Hence we would have two expla- 
nations for this shift: (a) the fact that regions with low 
strain (e.g vorticies) have smaller velocities than regions 
of high strain; (b) the existence of barriers of transport 
at the edge of the vorticies that act like parcel traps. 

We have shown that the distribution of r can be ob- 
tained from the strain distribution below one turnover 
time. Its time evolution at longer time scale has to be 
studied in more detail, especially its asymptotic form. 
The dependence between r and A is a complex issue 
which needs to be addressed. We expect it to be strongly 
dependent on the nature of the flow, especially on its 
Lagrangian correlation time. In our two-dimensional 
Navier-Stokes flow, areas of low stretching associated 
with long correlation times (e.g. vorticies) could explain 
the strong dependence between A and - where they are 
both small compared to their ensemble mean (figure [1} . 

Some studies (e.g. [l^, [2^) have shown, both nu- 
merically and theoretically, the relevance of Lagrangian 
stretching theories when the length scale of the tracer 
are much smaller than the scale of the flow, like in the 
present study. However, this was done for simple pre- 
scribed smooth chaotic flows and in the long time decay. 
This theory was recently applied by [2l| to infinite chem- 
istry in the long time decay. Here we apply it for the 
initial regime. A comprehensive picture would be given 
by studying intermediate time scales, which will be the 
subject of a future paper. 
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